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Abstract 

Background 

Recent genome-wide association studies (GWAS) have identihed a large number of genetic risk factors for schizophrenia 
(SCZ) featuring ion channels and calcium transporters. For some of these risk factors, independent prior investigations 
have examined the effects of genetic alterations on the cellular electrical excitability and calcium homeostasis. In the 
present proof-of-concept study, we harnessed these experimental results for modeling of computational properties on 
layer V cortical pyramidal cell and identify possible common alterations in behavior across SCZ-related genes. 
Methods 

We applied a biophysically detailed multi-compartmental model to study the excitability of a layer V pyramidal cell. 
We reviewed the literature on functional genomics for variants of genes associated with SCZ, and used changes in 
neuron model parameters to represent the effects of these variants. 

Results 

We present and apply a framework for examining the effects of subtle single nucleotide polymorphisms in ion channel 
and Ca^“'" transporter-encoding genes on neuron excitability. Our analysis indicates that most of the considered SCZ- 
related genetic variants affect the spiking behavior and intracellular calcium dynamics resulting from summation of 
inputs across the dendritic tree. 

Conclusions 

Our results suggest that alteration in the ability of a single neuron to integrate the inputs and scale its excitability may 
constitute a fundamental mechanistic contributor to mental disease, alongside with the previously proposed deficits 
in synaptic communication and network behavior. 


Introduction 


Schizophrenia (SCZ) is a severe mental disorder with heritability estimates ranging from 0.6 to 0.8 [Ripke et ah, 201^ . 
A recent genome-wide association study (GWAS) has identified more than a hundred genes exceeding genome-wide 
significance, confirming the polygenic nature of this psychiatric disorder [Ripke et ah, 2014] . This remarkable success 
in gene discovery brings up the next big challenge for psychiatric genetics - translation of the genetic associations into 
biological insights [van Os and Kapur, 2009| . Attaining this goal is supported by the development of biophysically 
detailed neuron models, boosted by the recent launch of mega-scale neuroscience projects jCrillner, 2014| . These 
models make it possible to investigate SCZ disease mechanisms by computational means, ultimately aiming towards 
achieving better clinical treatments and disorder outcomes [Wang and Krystal, 2014 [Oweip_2014J^ 


The 108 recently confirmed SCZ-linked loci span a wide set of protein-coding genes Ripke et ah, 2014], including 


numerous ion channel-encoding genes. The disorder is associated with genes affecting transmembrane currents of all 
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major ionic species, Na'*', K+, and Ca^’*'. In addition, some of the SCZ-linked genes are involved in regulation of the 
Cd?'^ concentration in the intracellular medium [Ripke et ah, 2014 , which is another great contributor to excitability. 
It is thus reasonable to hypothesize that the SCZ-linked genes should have an impact on the excitability at the single 
neuron level. 

We focus our study on cortical layer V pyramidal cells (L5PC) as a principal computational element of the cerebral 
circuit. An L5PC extends throughout the cortical depth with the soma located in layer V and the apical dendrite 
branching into the “apical tuft” in layer I, and its long axon may project to non-local cortical and subcortical areas. 
The tuft serves as an integration hub for long-distance synaptic inputs, and is often considered a biological substrate for 
cortical associations providing high-level “context” for low-level (e.g., sensory) inputs to the perisomatic compartment 
[Larkum, 20I3| . Therefore, the ability of L5PC to communicate the apical inputs to the soma has been proposed as one 
of the mechanisms that could be impaired in the mental disease [Larkum, 20I3| . In agreement with this hypothesis, 
recent psychiatric GWASs consistently reported association of genes coding for the subunits of voltage-gated Ca^’*' 
channels as risk factors in SCZ and bipolar disorder [Sklar et ah, 2011 Ripke et ah, 201 ij Smoller et ah, 2013] . 

In the present proof-of-principle study, we apply a model Hay et ah, 2011| of L5PC to explore how genetic variants 
in SCZ-linked genes affect the single-cell excitability. We carry out our study by linking a documented effect of a 
genetic variant in an ion channel or Ca^+ transporter-encoding gene to a change in the corresponding neuron model 
parameter. It should, however, be noted that information does not generally exist for the effect of single nucleotide 
polymorphism (SNP) variants identified through GWASs on the biophysical parameters required for the computational 
models. We instead use information obtained from in vitro studies of more extreme genetic variations, including loss 
of function mutations. A central assumption of this approach is that the effects of SNP variants can be represented 
as scaled-down versions of those of the more extreme variants, and that the emergence of the full psychiatric disease 
phenotype results from the combined effect of a large number of subtle SNP effects [Gottesman and Shields, 1967 


Purcell et ah, 2009 


Fromer et al., 2014 


A deficit in synaptic communication is likely to contribute to SGZ Stefansson et al., 2009 


Wen et ah, 2014 , Ripke et al., 2014 but is outside the scope of the present work. 


Methods and materials 

The L5PC model 

The multi-compartmental neuron model used in this work is based on a reconstructed morphology of a layer V thick- 
tufted pyramidal neuron (cell #1 in [Hay et ah, 201 1| ). The model includes the following ionic currents: Fast inactivat¬ 
ing Na’*' current (iNat), persistent Na"*" current (iNap), non-specific cation current (Ih), muscarinic K+ current (Im), 
slow inactivating K’*' current (Ikp), fast inactivating K"*" current {iKt)^ fast non-inactivating K+ current {Ikv 3 .i), high- 
voltage-activated Ga^"*" current (IcaHVA), low-voltage-activated Ca^''" current (IcaLVA), small-conductance Ca^+- 
activated K"*" current {IsK)^ and finally, the passive leak current (Iieak)- See Supplemental information for the model 
equations, and for the simulation codes, see ModelDB entry 169457 (https://senselab.med.yale.edu/ModelDB). 

Genes included in the study 

Since we did not aim to provide a comprehensive evaluation of a representative fraction of the genetic risk factors of 
SGZ, but to provide a proof of principle of the computational modeling approach, we selected the genomic loci using 
the following approach. We based our study on the recent GWAS [Ripke et al., 2014] , which reported significantly 
associated SNPs that were scattered across hundreds of genes with a variety of cellular functions. We concentrated on 
those genes that encoded either ionic channels or proteins contributing to transportation of intracellular Ca^“'" ions. 

We used the SNP-wise p-value data of [Ripke et al., 2014] , and for each gene of interest, determined the minimum 
p-value among those SNPs that were located in the considered gene. We performed this operation for all genes encoding 
either subunits of voltage-gated Ca^+, K+, or Na"*" channel, subunits of an SK, leak, or hyperpolarization-activated 
cyclic nucleotide-gated (HGN) channel, or Ca^+-transporting ATPases. The genes CACNAIC, CACNB2, CACNAlI, 
ATP2A2, and HCNl possessed a small minimum p-value each (p < 3 x 10“®) — these genes were also highlighted in 
the locus-oriented association analysis as performed in [Ripke et ah, 2014] . In order to extend our study to explore 
a larger set of genes, we used a more relaxed threshold (p < 3 x 10“®) for the minimum p-value, and obtained the 
following genes in addition to the previously mentioned ones: CACNAlD, CACNAlS, SCNIA, SCN7A, SCN9A, 
KCNN3, KCNS3, KCNBl, KCNG2, KCNH7, and ATP2B2. Of these, we omitted the genes that are not relevant 
for the firing behavior of an L5PG. It should be noted that we used the SNPs reported in [Ripke et ah, 2014] only 
to name the above SGZ-related genes, and due to lack of functional genomics data, we could not include the actual 
SCZ-related SNPs in our simulation study. Instead, we searched in PubMed for functional genomic studies reporting 
the effects of any genetic variant of the above genes. For details, see Supplemental information. 
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Results 


A new framework for bridging the gap between GWASs of SCZ and computational 
neuroscience 


In this work, we reviewed the literature on effects of variants in SCZ-related genes on ion channel behavior and 
intracellular Ca^“'" dynamics, and interpreted the reported effects in the context of our neuron model parameters. An 
overview of the relevant studies is given in Table while the effects of each variant on the L5PC model parameters 
are given in Table These data gave us a direct interface for linking a change in the genomic data, such as a SNP or 
an alternative splicing, into a change in neuron dynamics. The reported data, however, often corresponded to variants 
with large phenotypic consequences that in general are absent in SCZ patients. To simulate subtle cellular effects 
caused by the common SNPs related to SCZ [Lee et ah, 2012 , we downscaled the variants of Tableby bringing the 
changed parameters closer to the control neuron values when the reported change caused too large an effect in the 
neuron firing behavior. Our approach is illustrated in Figure 


Table 1: Table of the genetic variants used in this study. For more details, see 

Tabled 

Gene Refs. 




Type of variant 

Cell type 

CACNAlC 

Kudrnac et al., 2009 


L429T, L434T, S435T, S435A, S435P 

TSA201 

CACNAlC 

Kudrnac et al., 2009 


L779T, I781T, I781P 

TSA201 

CACNB2 

Cordeiro et al., 2009 


T Ill-mutant 

TSA201 

CACNB2 

Massa et al., 1995 



A1B2 vs Al alone 

HEK293 

CACNB2 

Link et al., 2009 



N1 vs N3 vs N4 vs N5 

HEK293 

CACNAID 

Tan et al., 2011 

•) 

Bock et ah, 2011 

42A splices transfected 

TSA201/HEK293 

CACNAID 

Tan et al., 2011 

•) 

Bock et ah, 2011 

43S splices transfected 

TSA201/HEK293 

CACNAID 

Zhang et al., 2011 

? 


Homozyg. knockout 

AV node cells / 


Perez-Alvarez et al., 2011| 


chromaffin cells 

CACNAII 

Murbartian et al., 2004 

Alternative splicing of exons 9 and 33 

HEK293 

ATP2A2 

Ji et al., 20001 




Heter. null mutation 

myocytes 

ATP2B2 

Fakira et al., 2012 



Heter. knockout 

Purkinje cells 

ATP2B2 

Empson et al., 2010 


Homozyg. knockout 

Purkinje cells 

ATP2B2 

Ficarella et al.. 

20071 

G283S-,G293S-mutant 

CHO cells 

SCNIA 

Cestele et al., 2008 


FHM mutation Q1489K 

Cultured neocortical cell 

SCNIA 

Vanmolkot et al., 2007 

FHM mutation L1649Q 

TSA201 

SCN9A 

Estacion et al.. 

2011 


I228M NaVl.7 variant 

HEK293 

SCN9A 

Estacion et al., 2008 


A1632E NaVl.7 mutation 

HEK293 

SCN9A 

Han et al., 2006 




L858F NaVl.7 mutation 

HEK293 

SCN9A 

Dib-Hajj et al., 2005 


F1449V NaVl.7 mutation 

HEK293 

KCNS3 

Shepard and Rae, 1999 

hKv2.1-(G4S)3-hKv9 fusion inserted 

HEK293 

KCNBl 

Bocksteins et al^^^^T 

T203K, T203D, S347K, S347D, T203W, S347W 

LTK- 

KCNN3 

Wittekindt et al^^OM 

hSK3_ex4 isoform 

TSA 

HCNl 

Ishii et ah, 2007 



D135W, D135H, D135N mutants 

HEK293 


The variants in Table were 23 in total, although some of them represented a range of effects of different variants 
(see Supplemental information. Table SI). The entries corresponded to variants of genes encoding for Ca^"*" chan¬ 
nel subunits (CACNAlC, CACNB2, CACNAID^ CACNAlI), intracellular Ca^’*' pumps {ATP2A2, ATP2B2), Na'*' 
channel subunits {SCNIA, SCN9A), K+ channel subunits {KCNS3, KCNBl, KCNN3), and a non-specific ion channel 
subunit (HCNl). In the following, we present simulation results for the L5PC model equipped with some of these 
downscaled variants (we refer to these model neurons as “variant neurons” or simply as “variants”). As we do not 
know the quality of the effects of the actual SCZ-related polymorphisms, we perform the simulations for a range of 
differently scaled variants, including negative scalings (i.e. opposite effects w.r.t. the effects reported in Table [^. We 
concentrated our study on a representative sample of six variants, highlighted in TableThese variants represent six 
genes with different roles in L5PC electrogenesis and a wide range of observed effects (see Table S2). 

Variants show altered intracellular [Ca^+] responses to short stimuli 

To characterize the implications of SCZ-related genes on the neuron excitability, we started by analyzing the effects 
of the downscaled versions of genetic variants in Table on the neuron response to a short somatic suprathreshold 
square-pulse stimulus. Figure IK -B shows the time course of the membrane potential and the intracellular Ca^^- 
concentration ([Ca^+J) for one variant, whereas Figure]^ shows the [Ca^+] response in the phase plane for several 
different variants. The most typical effect of a variant was a deviation in the peak [Ca^+], but differences could also 
be observed in the rising phase of the [Ca^"*"], as shown in the phase plane representation in Figure]^. 

The results in Figure]^ show that the [Ca^+J response served as a more sensitive indicator of genetic effects than 
the membrane potential. In CACNAlC and CACNB2 variants, which affect the high-voltage-activated (HVA) Ca^"*" 
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Figure 1: An example of downscaling with a variant [Link et al., 200^ of the gene CACNB2. The authors 
of [Link et ah, 2009] transfected different variants of CACNB2 DNA into HEK293 cells together with CaVl.2 subunit 
DNA, and measured the activation and inactivation curves of the Ca^'*' currents through the cell membrane. The 
cells transfected with different variants showed different values of mid-points and time constants for channel activation 
and inactivation: The mid-points of activation and inactivation varied by ±4.9 and ±5.1 mV among the variants, 
respectively, while the time constants of activation and inactivation varied from —40% to ± 68 % and —40% to ± 66 %, 
respectively. Here we illustrate a variant representing one possible combination of the endpoints of these ranges. 
A: The voltage-dependence of steady-state activation (moo) and inactivation (hoc) according to the neuron model 
|Hay et ah, 201 1| . Red curves represent the unsealed variant: Activation parameters VoSma and 14ffmb are changed 
by Al = —4.9 mV, inactivation parameters Voffha and Voffhb are changed by A 2 = 5.1 mV, activation time constants 
Tina and Tmb are increased by 68 % (773 = 1 . 68 ), and inactivation time constants T^a and Thb are increased by ± 66 % 
(774 = 1.66). Purple curves show the downscaled e = 5 variant (see below), and blue curves show the control neuron 
activation properties. B— C: The voltage-dependence of time constants for activation (B) and inactivation (C). D— J: 
Illustration of the scaling conditions. The downscaling is based on five conditions (see Supplemental information), 
none of which should be violated in the downscaled variant. In panels (D)-(F) (conditions I-III), the variant neuron 
should respond with the same number of spikes as the control neuron (blue). In panel (I), the integrated difference 
between the variant and control f-I curves should not exceed 10% of the integral of the control curve (condition IV). 
In panel (J), the difference between the limit cycles should not exceed a set limit (condition V). It can be observed 
that the unsealed variant (red) violates the conditions I and IV. To scale down the variant, each applied parameter 
change is brought to a fraction c < 1 , all in proportion, until the threshold of violating/non-violating variant is found. 
In this example, the threshold variant corresponds to parameter changes c x Ai, c x A 2 , 773 , and 77 I, where c = 0.4574. 
Any non-negative value below the threshold value c yields a downscaled variant that obeys conditions I-V. The purple 
curves represent the variant corresponding to parameter changes ce x Ai, ce x A 2 , 77 “, and 77 “, where c = 0.4574 and 
e = ^. The blue curves show the control neuron firing behavior. (D): Somatic membrane potential as a response to a 
150ms somatic square-pulse stimulus (0.696 nA). (E): Somatic response to a distal apical synaptic conductance with 
maximum 0.0612 qS. (F): Somatic response to a combination of somatic square-pulse current (10 ms x 1.137 nA) 
and a distal apical synaptic conductance (maximum 0.100 qS). (G): Somatic membrane potential as a response to a 
long somatic square-pulse stimulus current (1.0 nA). (H): Somatic Ca^"*" concentration as a response to the stimulus 
used in (G). (I): Spiking frequency as a function of amplitude of the somatic DC. (J): The membrane potential limit 
cycle corresponding to the late phase of (G). The x-axis represents the somatic membrane potential and y-axis its 
time derivative. a 



























































Table 2: List of variants and their threshold effect coefficients c. The variants are ordered as in Table but 
the variants where several combinations of parameter range end points were considered are listed as separate variants 
(see Supplemental information). The variants marked with f were used in Figures |2j|^ the variants marked with § 
were used in Figure S3, and the variant marked with | was used in Figure [T| 


Gene 

Effect on model parameters 

Thresh, effect 

CACNAIC 

Voffm*,CaHVA: -25.9 mV; Voffh*,CaHVA: -27.0 mV 

c ^ 0.094 

CACNAIC 

Voffm*,CaHVA: -37.3 mV; Voffh*,CaHVA: -30.0 mV 

c = 0.060 t 

CACNB2 

Voffh*,CaHVA: -5.2 mV; Vsioh*,CaHVA: xO.69 

c ^ 0.582 

CACNB2 

Th*. CaHVA: Xl.7 

c > 2.000 


CAGNB2 


^^ffm* ,CaHVA 
V'offm*,CaHVA 
^offm-x.CaHVA 
V'offm*,CaHVA 
V'offm-x.CaHVA 
Voffm*,CaHVA 
V'offm*,CaHVA 
^offm*,CaHVA 
Voffm*,CaHVA 
V'offm-x.CaHVA 
V'offm*,CaHVA 
^offm*,CaHVA 
^offmK.CaHVA 
V'offm*,CaHVA 
V'offm-x.CaHVA 
^offm-x.CaHVA 


-4.9 mV; Voffh*,CaHVA: -5.1 mV; Tm^.CaHVA: xO.6; rh*,CaHVA: XO.6 
-1-4.9 mV; Voffh*,CaHVA: -5.1 mV; Tm^^^CaHVA: xO.6; Th*,CaHVA: xO.6 
-4.9 mV; Voffh*,CaHVA: -h5.1 mV; Tin.^,CaHVA: xO.6; Th^,CaHVA: xO.6 
-1-4.9 mV; Voffh*,CaHVA: -1-5.1 mV; Tm*,CaHVA: xO.6; Th*,CaHVA: xO.6 
-4.9 mV; Voffh*,CaHVA: -5.1 mV; rni*,CaHVAi xl.68; ThK.CaHVA: xO.6 
-1-4.9 mV; Voffh*,CaHVA: -5.1 mV; Tm*,CaHVA: Xl.68; Th..,CaHVA: xO.6 
-4.9 mV; Voffh*,CaHVA: -K5.1 mV; rin»,CaHVA: xl.68; Th.s,CaHVA- xO.6 
-1-4.9 mV; V)ffh-^,CaHVA: -1-5.1 mV; Tm*,CaHVA: Xl.68; Th»,CaHVA: xO.6 
-4.9 mV; Voffh*,CaHVA: -5.1 mV; rm*,CaHVA: xO.6; rh*,CaHVA: Xl.66 
-1-4.9 mV; Voffh*,CaHVA: -5.1 mV; rxn*,CaHVA: xO.6; Th*,CaHVA: xl.66 
-4.9 mV; V'offh*,CaHVA: -h5.1 mV; Tm*,CaHVA: xO.6; Th*,CaHVA: Xl.66 
-1-4.9 mV; V^ffh*,CaHVA: 4-5.1 mV; Tm*,CaHVA: xO.6; T’h*,CaHVA: xl.66 
-4.9 mV; Voffh*,CaHVA: -5.1 mV; Tni*,CaHVA: xl.68; Th+,CaHVA: Xl.66 
-1-4.9 mV; Voffh*,CaHVA: -5.1 mV; Tm*,CaHVA: Xl.68; Th..,CaHVA: xl.66 
-4.9 mV; VoffhH<,CaHVA: -45.1 mV; rxn*,CaHVA: xl.68; Th«,CaHVA: xl.66 
-1-4.9 mV; V^ffh*,CaHVA: 4-5.1 mV; Tm*,CaHVA: Xl.68; Th*,CaHVA: xl.66 


c = 0.310 
c ^ 0.154 
c ^ 0.205 
c ^ 0.675 § 

c ^ 0.814 
c ^ 0.057 
c ^ 0.517 
c ^ 0.078 
c ^ 0.275 
c ^ 0.188 
c ^ 0.190 
c ^ 1.687 
c ^ 0.707 
c ^ 0.061 
c ^ 0.457 t t 
c = 0.084 


CACNAID 


VoffmK 

K>ffmK 


CaHVA 

,CaHVA 


-10.9 mV; 
-10.9 mV; 


Klom 

^lom 

Klom 

14lom 

V.lom 

^lom 

^slom 

l^lom 


CaHVA 

CaHVA 


xO 

xO 


■ 73; Voffh*,CaHVA: -3.0 mV; Vsioh*,CaHVA: x0.81; Th*,CaHVA: xl.25 
■73; Voffh*,CaHVA- 4-3.5 mV; Vsioh*,CaHVA: x0.81; Th*,CaHVA: xl.25 


c = 0. 
c = 0. 


118 

106 


CACNAID 


Voffm-K 

^^ffm* 

^^ffm* 

Voffm* 

Voffm* 

^^ffm* 

^^ffm* 

V>ffm* 


CaHVA 

CaHVA 

CaHVA 

CaHVA 

CaHVA 

CaHVA 

CaHVA 

CaHVA 


-10.6 mV: 
-43.4 mV; 
-10.6 mV: 
-43.4 mV; 
-10.6 mV: 
-43.4 mV; 
-10.6 mV: 
-43.4 mV; 


CaHVA 
CaHVA 
*, CaHVA 
CaHVA 
CaHVA 
*, CaHVA 
=tc, CaHVA 
*, CaHVA 


xO.8; Voffh*, CaHVA: -5.3 mV; Vsioh*, CaHVA: 
xO.8; Voffh*,CaHVA: -5.3 mV; Vsioh*,CaHVA: 
xl.12; Voffh*,CaHVA: -5.3 mV; V^loh*,CaHVA 
xl.12; Voffh*, CaHVA: -5.3 mV; Vsloh*, CaHVA 
xO.8; Voffh*, CaHVA: 4-1.2 mV; Vsioh*, CaHVA 
xO.8; Voffh*,CaHVA: -41-2 mV; Vsioh*,CaHVA 
xl.12; Voffh*,CaHVA: 4-1.2 mV; Vsioh*,CaHV 
Xl.12; Voffh*. CaHVA: 4-1.2 mV; Vsioh*. CaHV. 


XO.66; Th*,CaHVA: 

xO.72 

c = 0.114 

xO.66; Th*,CaHVA: 

xO.72 

c ^ 1.962 

XO.66; Th*,CaHVA 

xO.72 

c ^ 0.131 

XO.66; Th*,CaHVA 

xO.72 

c ^ 0.601 

XO.66; Th*,CaHVA 

xO.72 

c ^ 0.100 

XO.66; Th*,CaHVA 

xO.72 

c ^ 0.645 

j.: xO.66; Th*,CaHVA 

: xO.72 

c ^ 0.116 

: xO.66; Th*,CaHVA 

: xO.72 

c = 1.117 


CACNAID 


V>ffm* 

Voffm* 

Voffm* 

Voffm* 


CaHVA 

CaHVA 

CaHVA 

CaHVA 


4-6.6 mV; 
4-6.6 mV; 
4-6.6 mV; 
4-6.6 mV; 


^lom* 

Vslom* 

14lom* 

14lom* 


CaHVA 

CaHVA 

CaHVA 

CaHVA 


xO. 

xl. 

xO. 

xl. 


75; Th*,CaHVA 
19; Th*,CaHVA 
75; Th*,CaHVA 
19; Th*,CaHVA 


xO.5 

xO.5 

xl.12 

xl.12 


c = 0. 
c ^ 0. 
c ^ 0. 
c = 0. 


104 

068 

115 

072 


CACNAII 


V>ffma, 

V>ffma, 


CaLVA: 
CaLVA: 


-41.3 mV; 
-41.3 mV; 


Vaffha. 

Voffha, 


CaLVA: 

CaLVA: 


+ 1.6 
+ 1.6 


mV; Tm*,CaLVA: XO.87; Th*,CaLVA: XO.i 
mV; Tm*,CaLVA: Xl.45; Th„^CaLVA: xO.; 


c > 2, 
c > 2. 


000 

000 t 


ATP2A2 _7CaDynamics: xO.6_c — 0.093 § 


ATP2B2 

Tdecay, CaDynamics • Xl.97 

c = 0.218 

ATP2B2 

Tdecay, CaDynamics • Xl.5, Cmin,CaDynamics • Xl.4 

c = 0.215 t 

ATP2B2 

Tdecay, CaDynamics • X4.45 

c ^ 0.099 

SCNIA 

Voffm,Nat: -0.3 mV; Vjffh,Nat: +5 mV; V^iom,Nat: Xl.l5; V^ioh,Nat: Xl.23 

c = 0.056 t 

SCNIA 

Voffm,Nat: +2.8 mV; Voffh,Nat: +9.6 mV; Vsiom,Nat: xo.984; VBioh,Nat: X 1.042 

c ^ 0.069 

SCN9A 

Voffh*,Nap: +6.8 mV 

c > 2.000 

SCN9A 

Voffh*,Nap: +3.5 mV; Vsloh,Nap: XO.55; Voffm,Nat: -7.1 mV; Voffh,Nat: +17.0 mV; Vloh,Nat: xO.69 

c ^ 0.026 

SCN9A 

Voffm,Nat: -9.1 mV; Voffh,Nat: +3.1 mV 

c = 0.043 

SCN9A 

Voffm,Nat: -7.6 mV; Voffh,Nat: +4.3 mV 

c = 0.043 

KCNS3 

Tm*,Kp: x2.0; Th*,Kp: X2.5; Vioh,Kp: xO.5 

c > 2.000 

KCNBl 

Voffm,Kp: +5 mV; Voffh,Kp: +3 mV; l+iom,Kp: xl.ll; Vloh,Kp: xO.86; Tm*,Kp: xO.5; Th*,Kp: xO.53 

c > 2.000 

KCNBl 

Voffm,Kp: +1 mV; Voffh,Kp: -0 mV; V^iom,Kp: xl.22; l+ioh.Kp: xl.O; Tj^,*,j^p: xO.89; Th*,Kp: Xl.l3 

c > 2.000 

KCNBl 

Voffm,Kp: +6 mV; Voffh,Kp: -8 mV; Viom.Kp: Xl.33; Vgioh.Kp: Xl.O; Tm*,Kp: xO.5; Th*,Kp: xO.87 

c > 2.000 

KCNBl 

Voffm,Kp: -28 mV; Voffh,Kp: -27 mV; Vgiom.Kp: xl.ll; Vsioh,Kp: x0.71; Tm*,Kp: Xl.l3; Th*,Kp: x2.27 

c > 2.000 

KCNBl 

Voffm,Kp: +14 mV; Voffh,Kp: -21 mV; Vsiom.Kp: x2.0; Vsioh,Kp: Xl.O; Tm*,Kp: xO.39; Th*,Kp: Xl.2 

c > 2.000 

KCNBl 

Voffm,Kp: -13 mV; Voffh,Kp: -13 mV; VBiom,Kp: xl.33; V^ioh,Kp: x0.71; Tm*,Kp: xO.95; Th*,Kp: x5.13 

c > 2.000 

KCNN3 

Coff,SK: xO.86; Cs 1 o,sk: Xl.24 

c = 1.715 t § 

HCNl 

Voffma,Ih: Voffmb,Ih: -26.5 mV; l+loma,Ih) l^lomb,Ih: XO.64 

c ^ 0.296 

current, and in CACNAII variant affecting the low-voltage-activated (LVA) Gs?^ current, the observed effects on 
the [Ca^+J response were due to changed Ca^’*' channel kinetics. Similarly, in the ATP2B2 variant, which affects 
the plasma membrane Ca^’*' ATPase (PMCA) activity, the observed effect was caused by altered intracellular Ca^’*' 
dynamics. By contrast, the small yet observable differences between control and SCNIA variant [Ca^’*'] phase planes 


were due to alterations in subthreshold membrane potential fluctuations, which caused variation in the activation of 
Ca^’*' channels. 


Steady-state firing is influenced by the variants 

Next, we investigated the steady-state behavior of the neuron when a direct current (DC) was applied to the soma. 
As shown in Figure]^ the f-I curves (firing frequency as a function of DC amplitude) of many SCZ-associated variants 
were notably different from those obtained with the control neuron. 

The deviations in the f-1 curves in many of the variants can be explained by changes in the Ca^“'"-activated SK 
current. These changes were caused either by direct alteration of the activation kinetics of SK channels (as in the 
KCNN3 variant), or indirectly, through the changes in the intracellular [Ca^+] response associated with other variants 
{CACNAIC, CACNB2 and ATP2B2). As an example, the Ca^+ influx during an action potential was larger in the 
CACNAIC variant neuron compared to the control neuron (compare curves in blue and magenta in Figure]^), and 
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ATP2B2 SCN1A KCNN3 
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Figure 2: Variants show altered [Ca^+] response to a short, somatic stimulus. A,B: The membrane potential 
(A) and Ca^"*" concentration (B) time courses, recorded at soma, as a response to a somatic 5 ms, 1.626 nA square- 
pulse stimulus. The blue curve shows the control neuron behavior, while the other curves show the behavior of 
a CACNAlI variant [Murbartian et ah, 2004 with different scalings (magenta: e = 5 , dark gray: e = |, cyan: 
e = — j, green: e = — ^). Inset: Zoomed-in view around the time of the spike. No notable differences between the 
variants and the control can be observed in the membrane potential time course. C: The [Ca^"*"] response plotted 
in the phase plane for different variants. The behavior of variants of CACNAlC [Kudrnac et ah, 2009 , CACNB2 


[Link et al., 2009 , CACNAlI Murbartian et ah, 2004| , ATP2B2 [Empson et ah, 2010| , SCNlA [Cestele et ah, 20081 , 
and KCNN3 Wittekindt et al., 2004 genes are shown with similar scaling as in (A) and (B). 


this led to an increase in the SK current activation, relative to that in the control neuron. The increased Isk, in 
turn, delayed the induction of the next spike, and hence resulted in a loss of gain (flattening) in the f-1 curve (Figure 
[^. This finding is in line with previous modeling studies, which also have highlighted the role of Ca^+-activated K+ 
currents in modulating f-I curves [Engel et al., 1999] Haines et ah, 2011]. 


— control — 6 = 1/2 —6 = 1/4 — 6 = -1/4 — 6 = -1/2 


CACNA1C CACNB2 CACNA11 



0.4 0.8 1.2 0.4 0.8 1.2 0.4 0.8 1.2 


I (nA) I (nA) I (nA) 


Figure 3: The variant neurons show modulated gain. f-I curves are shown for different scalings of different 


variants (magenta: e = 5 , dark gray: e = j, cyan: e = — |, 
CACNAlC, CACNB2, ATP2B2, and KCNN3 variants. 


green: e = —h)- Differences in gain are visible for 
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We also studied the steady-state firing behavior by analyzing the limit cycle, i.e., the phase plane representation 
of the variable of interest after a large number of initial cycles. As the [Ca^’'"] response served as the most sensitive 
biomarker, we used the [Ca^’*'] limit cycle for this analysis. Figure]^ shows, for some of the gene variants, the [Ca^’*'] 
limit cycle as a response to a long, square-pulse stimulus. 



Figure 4: Differences in [Ca^+] limit cycles of the variants. A,B: The membrane potential (A) and Ca^'*' 
concentration (B) time courses at soma as a response to a DC with an amplitude of 1.2 nA. The different colors 
represent different scalings of a CACNAlC variant [Kudrnac et ah, 2009| (scaling as in Figure]^. Inset: The time 
course shown 3.4 s since the beginning of the stimulus. C: The intracellular [Ca^’*'] phase plane during steady firing 
caused by a DC applied to the soma. The variants and their scalings were chosen as in Figure]^. 

The different variants had different ways of impacting the shape of the [Ca^"*"] limit cycle. An especially common 
trait was a shift or compression in the horizontal direction, representing a small increase or decrease of intracellular 
[Ca^^] during steady hring. Figure SI shows the range of [Ca^"*"] values observed during the constant current injection 
for all variants used in this study. Interestingly, almost all variants gave deviations from the [Ca^"*"] range observed in 
the control neuron. 

The only genes that did not show notable variance for the range of [Ca^"*"] during steady-state hring were CACNAlI, 
KCNS3, and KCNBl. Despite this, CACNAlI did have a clear effect on the [Ca^+J phase plane in the single-pulse 
response (see Figure [^). Conversely, the KCNS3 and KCNBl variants were not found to inhuence any of the 
considered aspects of neuron behavior (see Figures SI and S2). In these variants, the activation and inactivation of 
the persistent K"*' current {Ikp) have been modihed. Our hndings rehect the overall minor role that the current Ixp 
plays in shaping the response properties of the L5PC model neuron: It is only present in the soma and with relatively 
small maximal conductance [Hay et ah, 2011] . The minor contribution of this current to the model neuron behavior 
may be in contradiction with experimental evidence, which states that this current constitutes a major proportion of 
the total outward K+ currents |Guan et ah, 2007] . 

Variants have an effect on integration of somatic and apical inputs 

A question that arises from the observation of the above trends is whether, and to what extent, the small deviations 
from control neuron behavior affect the information processing capabilities of the neuron. In [Hay et ah, 20 ll] , it was 
shown that the L5PC model can be used to describe the Ca^’*' spike generation as a response to a combination of 
stimuli at the apical dendrite and at the soma. Moreover, they showed a good qualitative match with experimental 
data and model predictions for the neuron responses during an “up” or “down” state. In this work, we followed their 
dehnition for the down state as the resting state of the neuron, and the up state as a state where a subthreshold 
current is applied to the proximal apical dendrite |Hay et ah, 201 1| . We also adopted their protocol of combining an 
apical and somatic stimulus and studying the effect of order and inter-stimulus interval (ISI) between the two (see 
Figure 8 in [Hay et ah, 201 1| ). This gave us a temporal profile showing a range of ISI for which a large Ca^’*' spike is 
produced. We then estimated the effect of our genetic variants on this temporal profile. 

Figure shows the temporal profiles for the six representative variants. Effects of the different variants on the 
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neuron response were visible both during the down and during the up state. Particularly large effects were observed in 
CACNAII and SCNIA variants, which code for the LVA Ca^’*' channels and transient Na"*" channels that contribute 
to the sharp rise of voltage during an action potential. These channel types are expressed both in the soma and in the 
apical dendrite. Although all downscaled variants shared the basic form of the response curves, the shape was clearly 
altered in some variants. As shown in [Hay et ah, 201 1| , the up-state temporal profile of the control neuron showed 
elevated apical responses when the apical stimulus was applied shortly after the somatic stimulus, but not when the 
apical stimulus was applied first. This order-specific response was changed in some of the variants (see e.g. SCNIA 
variant) so that they also produced large Ca^+ spikes when the apical stimulus was applied shortly before the somatic 
stimulus, and could hence alter the order sensitivity of the coincidence detection. 


CACNA1C CACNB2 CACNA11 




ISI (ms) ISI (ms) ISI (ms) 


Figure 5: Sensitivity of variant neurons to ISI between somatic and distal apical stimuli during up and 
down states. The magnitude of the Ca^+ spikes was assessed by measuring the simulated membrane potential near 
the bifurcation point of the apical dendrite, i.e, at a distance of 620 pm from the soma. The dashed lines show the 
temporal maximum of the model response membrane potential during a down state, and the solid lines during an up 
state. The x-axis represents the ISI between the somatic and apical stimuli, positive values denoting cases where the 
apical stimulus was applied after the somatic stimulus. The variants shown are the same as in Figure In the up 
state paradigm, the neuron was first given a depolarizing square-pulse current of 0.42 nA x 200 ms at the proximal 
apical dendrite, 200 pm from the soma. In the middle of this period, a somatic square-pulse current of 0.5 nA x 5 ms 
was applied, and after a time defined as the ISI, an alpha-shaped current injection with rise time 0.5 ms, decay time 
5 ms, and maximal amplitude of 0.5 nA was applied. The down state paradigm was otherwise equal, but the long 
depolarizing current was absent, and to compensate this, the short somatic square-pulse current had an amplitude of 
1.8 nA. 


Variants show differences in inhibition of a second apical stimulus 

The alteration of information processing by the variants can also be seen in the sensitivity to successive synaptic 
inputs at the same locations. To explore this, we gave two successive stimuli to the apical dendrite, and varied the 
interval (ISI) between them. The first stimulus activated long time scale inhibitory currents, especially Isk and Im, 
and within some following time window, these currents could inhibit the generation of a subsequent spike. This could 
be interpreted as a form of single-cell pre-pulse inhibition: If the initial input caused the neuron to spike, the next 
input could fail to do so, even at times when the second stimulus were of larger amplitude. 

Figure shows how the threshold synaptic conductance needed for inducing a second spike, relative to the one 
needed for the first spike, depends on the ISI. Interestingly, the threshold conductance curve was affected by many of 
the variants, especially by the variants of Ca^'*' channel genes. Although the shown effect was partly attributed to the 
differences in the spiking thresholds between the variants, differences between the variants could still be observed when 
an absolute spiking threshold values for inducing a second spike were considered (data not shown). These findings 
raise the possibility that the disturbed prepulse inhibition observed in many SCZ patients might be an effect of altered 
ionic channel properties at the single neuron level, and not only of the modified elements in synaptic circuitry as 
previously thought [Swerdlow and Geyer, 1993 Medan and Preuss, 2011| (cf. [Koch, 20121). 


Combinations of variants may cause large alterations of neuron excitability 

The effects of the variants on neuron excitability, as explored in Figures [2]|^ were relatively small when the variants 
were studied in isolation (this was also the intention behind the downscaling process). However, combinations of 
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Figure 6: Threshold conductance factor for inducing a second spike as a function of ISI. The neuron was 
first made to spike by applying 3000 simultaneously activated synaptic conductances that were uniformly distributed 
along the apical dendrite at a distance [300 pm, 1300 pm] from the soma. All synaptic conductances were alpha-shaped 
with r = 5 ms and maximal conductance l.lpthj where t/th was the threshold conductance value for inducing a spike in 
resting state. These values were variant-specific (and scaling-specific), and they are plotted in the inset of each panel. 
After the specified ISI, the synapses were activated again. The y-axis shows the factor Cg, chosen such that c^gth was 
the threshold conductance of the second stimulus for generating an extra spike. Curves above the control curve (blue) 
represent strengthened prepulse inhibition, while curves below it represent weakened prepulse inhibition. The effect 
of the KCNN3 variant (e = |, magenta) was bilateral: For ISI < 50 ms, the effect was strengthening, whereas for ISI 
> 50 ms, the effect was weakening. 


variants of different genes may have additive effects on the neuron excitability. As a final step of our inquiry, we went 
on to investigate the effect that combinations of several of the variants could have. Figure S3 shows a selected example 
of how a combination of variants altered the neuron firing properties considered in Figures [2]j^ Although this version 
of the model neuron consists of downscaled variants, the properties of neuron excitability were remarkably modified. 


Discussion 


We presented a framework for studying the effects of genetic variants in genes related to schizophrenia (SCZ) on layer 
V pyramidal cell (L5PC) excitability. Our results indicate that most of the studied variants predict an observable 
change in the neuron behavior (see Figures [2j^. Taken together, these data provide support for using neuronal 
modeling to study the functional implications of SCZ-related genes on neuronal excitability. Although the analyses 
presented in this paper are specific to L5PCs and SCZ-related genes, our framework may be directly applicable to 
other cell types and other polygenic diseases, such as bipolar disorder and autism, given an identification of risk genes 
related to neuronal excitability. Furthermore, our approach can be directly extended to biophysically detailed models 
of neuronal networks, e.g. |Hay and Segev, 20I4| . 

The results shown in the present paper are based on a multi-compartmental neuron model |Hay et ah, 201 1] em¬ 
ploying a certain set of ionic channels. Alternative models with comparable complexity have recently been published. 
The authors of |Bahl et ah, 2012] present and apply a method for fitting the L5PC model to a simplified neuron mor¬ 
phology. In another recent work |Almog and Korngreen, 2014 , the L5PC model was fit to each stained morphology 
and the corresponding recordings separately. To rule out our conclusions being an artefact of a certain specific morphol¬ 
ogy, we carried out our analyses on another neuron model fit (cell #I equipped with channel conductance values from 
the first column of Table SI in [Hay et ah, 201 1| ) and an alternative neuron morphology (cell #2 of [Hay et ah, 201 1| ). 
The trend of the variant effects remained the same in these simulations (see Supplemental information, Figures S4 
and S5). More comprehensive analysis could be carried out by using the above-mentioned alternative models, or by 
employing a wider set of neuron morphologies |Hay et ah, 2013] , but this is left for future work. 

Discovering the disease mechanisms of disorders with polygenic architecture, such as SCZ, requires considering 
the interplay between many different biological processes and entities. The computational framework presented here 
allows one to identify physiological mechanisms and biomarkers common across multiple risk-related genes. Although 
there are many more SCZ-related genes [Ripke et ah, 201^ than those listed in Table the identified variants capture 
multiple types of changes in neuron excitability (see Figures [2]-[^. This allows studying the effects of combinations of 
different genetic variants in a straightforward manner, as suggested in Supplemental information (Figure S3), although 
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one should be careful in the interpretation of the results as certain variants could have non-additive effects. The gene 
hits not considered in the present work contribute to e.g. expression level and function of ionotropic and G-protein 
coupled receptors, as well as phosphorylation and dephosphorylation of diverse proteins. It should be noted that 
several of the identified SCZ-related genes, such as PTN, PAK6, and SNAP91 may also affect neuron morphology 
[Deuel et ah, 2002][Nekrasova et ah, 2008|[Ang et ah, 2003[[Dickman et ah, 2006[ . Our framework could be extended 
to include some of these effects, but studying the contribution of many signaling genes is out of the scope of the used 
neuron model. 

Linking cellular mechanisms to SCZ can be extremely difficult due to the high cognitive level on which the disease 
symptoms usually emerge. There is however one important exception, namely, altered prepulse inhibition, which 
is found in approximately 40-50% of SCZ patients [Turetsky et ah, 2007[ . In the present study, we showed that 
the considered genetic variants may have different tendencies to suppress a second synaptic stimulus, owing only to 
the differences in the (non-synaptic) ion channel gating or Ca^^ dynamics (Figure]^. This could be an important 
contributor to understanding the phenomenon of prepulse inhibition, although deficiencies in synaptic connections 
are most likely to play a role as well [Swerdlow and Geyer, 1993) . Moreover, altered integration of local and distal 
inputs in the manner of Figure could possibly underlie incomplete or excessive activation of pyramidal cells in 
certain brain areas linked to auditory hallucinations ( [Larkum, 2013[ ; see also discussion of the role of spontaneous 
activity in primary auditory cortex in triggering verbal hallucinations in [Kompus et ah, 2013[ ). Experimental testing 
of the shown model predictions for integration of inputs on a single-cell level is possible as well. Even without genetic 
manipulations, the downscaled variants of the ion channels could be imitated in vitro by partial pharmacological blocks 
that are configured to cause changes in currents comparable to those in our variants. However, in such a case the fine 
details on e.g. current activation time constants would be dismissed. 

The genetic variation of ion channels will affect not only the membrane potentials of the neurons, but also 
the transmembrane currents generating the local field potential (LFP) recorded inside cortex and the electrocor- 
ticography (EGoG) signals recorded at the cortical surface [Buzsaki et ah, 2012 Einevoll et ah, 2013[ . Further, the 
same transmembrane currents also determine the neuronal current dipole moments measured in EEG and MEG 
[Nunez and Srinivasan, 2006[ [Hamalainen et ah, 1993[ . Thus a natural extension of the present work is to investi¬ 
gate the effects of the genetic variation of the single-neuron current dipole moments [Linden et ah, 20T0] and the 
corresponding single-neuron contribution to the EEG signals. This is left for future work. 

A severe limitation of our approach is the lack of data on how different specific SCZ-related genetic polymorphisms 
affect the ion channel function and Ca^“'" dynamics. In this paper, we used literature data on extreme variants (e.g., 
loss of function mutations) of the identified SCZ-related risk genes, from a range of cell types, including non-neural 
tissues. Ideally, one would want to use physiological measures for the relevant genetic variants in the actual cell types 
(i.e., L5PCs) in cortical circuits, but such data do not currently exist. On the other hand, the diversity of data in the 
literature that we used and the possible discrepancies in the experimental procedures therein could also reflect in the 
results shown in this study, although our downscaling procedure constrains such errors from above. Novel methods 
employing automated cell patching [Ranjan et ah, 2014] to test the effects of subtle SCZ-related SNPs could provide 
the solution to both above limitations of our study. 

The validity of the downscaling, which for the aforementioned reasons is a central procedure in our framework, is 
built upon two assumptions: 1) That the SCZ-related gene variants do affect neuronal excitability on a single-neuron 
level, and 2) that their effects on ion channel or calcium pump kinetics are smaller than those of certain already 
documented mutations with extreme phenotypic consequences, such as deafness, hemiplegic migraine, or other pain 
disorders. The first assumption can be argued for by the multitude of ion channel-encoding genes that have been shown 
to be related to SCZ risk, and the second by the lack of severe body malfunctioning (such as the abovementioned 
ones) in SCZ patients, but rigorous proofs for both assumptions are yet to be shown. A challenge for future studies 
is also to decipher how the ion channel densities are affected by the variants in various cell types, an issue that might 
be important in modeling the effects of SCZ-associated polymorphisms. 

Our framework is an early attempt toward understanding the disease mechanisms of polygenic psychiatric disorders 
by computational means. It binds together genetic and single neuron levels of abstraction in the big picture of modeling 
psychiatric illnesses, as recently sketched in [Wang and Krystal, 2014[ and [Corlett and Fletcher, 2014 . While earlier 


computational studies considered only effects of single variants on neuron firing behavior Murbartian et ah, 2004 


Zhang et ah, 20TT| , our framework allows the screening of multiple types of variants and their implications on the 
neuron excitability. It may serve as a proof of principle for a novel “Biophysical Psychiatry” approach, enabling 
the integration of information from previously disparate fields of study including GWASs, functional genomics, and 
biophysical computational modeling using models of moderate complexity. We consider the chosen level of complexity 
a suitable trade-off for this purpose: While more approximate models, such as integrate-and-fire models, could not 
distinguish between the effects of different genetic variants in an acceptable precision, more detailed models including 
e.g. dynamics of protein translation and phosphorylation are likely to be very computation-intensive. The challenge 
for future work is to extend the study to both larger and finer spatial and temporal scales. 
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Functional effects of schizophrenia-linked genetic variants on intrinsic 
single-neuron excitability: A modeling study. Supplemental information 


Methods 

The L5PC model 

The current balance equation for each compartment of the neuron can be written as 
dV 

= I Nat + I Nap + + Im + Ixp + ^Kt + IkvS.I + ICaHVA + ICaLVA + IsK + h + laxiah 

where each type of current, except for the axial current, can be described as a multiplication of activation and inactivation 
variables as 

I = gm^'^h^''{E -V). 

Here, g is the maximal conductance of the ion channels, m and h are the activation and inactivation variables while Nm and 
Nh are constants describing their sensitivities, and E is the reversal potential corresponding to the ionic species. Reversal 
potentials of Na+ and K+ are constants {Ef^a = 50 mV, E^ = —85 mV), while the reversal potential of Ca^’*' depends on 
the intracellular [Ca^+J: Eca varied between 96 mV and 120 mV in our simulations. The dynamics of the activation and 
inactivation variables are defined as 


dm m — moo dh h — hoo 

—— = - and — = -, 

dt Tm dt Th 

where moo, hoo, ^m, and Th are functions of membrane potential V. Typically, functions moo and hoo have a sigmoidal 
shape, where the half-activation and half-inactivation voltages are determined by one or more (depending on ion channel) 
parameters each. We denotate these parameters Vffm* and where * stands for further specifications if there are 

multiple parameters affecting the half-(in)activation voltage. In a similar fashion, parameters I4iom* and I4ioh* affect the 
slopes of the (in)activation curves, and parameters Tm» and Toffh* influence the time constants. As an exception, the activation 
of IsK is solely dependent on the intracellular [Ca^^], and this dependence is quantified by half-activation concentration 
parameter CoS and slope parameter Cgio. The intracellular [Ca^^] obeys the following dynamics: 

d[Ca^~*~]t _ IcaHVA + ICaLVA _ [Ca^~'~]i — Cmin 
dt 2.^E d rdecay 

where IcaHVA and IcaLVA are the high and low-voltage activated Ca^’*' currents entering the considered cell segment, 7 
represents the fraction of Ca^’*' ions entering the cell that contribute to the intracellular [Ca^+J, E the Faraday constant, 
d is the depth of the sub-membrane layer considered for calculation of concentration, Cmin the resting intracellular [Ca^'*'], 
and Tdecay Is the decay time constant of the intracellular [Ca^+J. The simulation codes are provided in the ModelDB entry 
169457 (https://senselab.med.yale.edu/ModelDB). 

Channel activation and inactivation dynamics 

Here, the model equations and parameter values are listed for all ion channels. The maximal conductances g are different 
for different parts of the neuron: The subindex s refers to soma, ad to apical dendrite, bd to basal dendrite, and ax to axon, 
and for each current type all unshown maximal conductances are zero. For the non-specific ionic currents and Ca^"*" currents, 
the maximal conductances vary spatially along the apical dendrite in the following way. In the “hot zone”, which lies on 
the distance 650|J.m-850|a.m, the maximal conductances of Ca^^ currents are ten or hundredfold larger than elsewhere in the 
apical dendrite. The maximal conductances of the non-specific ionic current in turn grow exponentially from 0 to the length 
of the longest branch, Tmax) which is 1300 qm in this study. The numerical values correspond to the control neuron values 
and may be varied in the variant neurons. 
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Non-specific cation current, Ih 
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Muscarinic K+ current, 
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High-voltage-activated Ca^+ current, IcaHVA 
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Low-voltage-activated Ca^+ current, IcaLVA 
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Functional genomics literature review 

Many of the genes that are now confirmed to be linked to SCZ, have previously been shown to play a role in regulating the 
activation/inactivation kinetics of certain types of cellular transmembrane currents in animal or in vitro studies. Typically, 
these studies involved transfection of ion channel-encoding DNA into cells that normally do not express the considered ion 
channels, and documentation on how a variant DNA changed the electrophysiological properties of the cells compared to the 
cells transfected with control DNA. Moreover, many studies have demonstrated the effects of certain variants of a calcium 
signaling toolkit gene [SI] on the dynamics of the Ca?'^ concentration in the intracellular medium. 

We searched through the literature on functional genomics for genes CACNAIC, CACNB2, CACNAlI, ATP2A2, HCNl, 
CACNAID, CACNAIS, SCNIA, SCN7A, SCN9A, KCNN3, KCNS3, KCNBl, KCNG2, KCNH7, and ATP2B2 to find data 
on how genetic variations change the ion channel behavior or intracellular Cd?'^ dynamics. Due to shortness of data reported 
for a single animal and cell type, we included studies performed in various animal species and across different tissues. We 
concentrated on studies that fulfilled the following conditions: 

• The study applied a genetic variant of one of the genes of interest and the gene is likely to be expressed in L5PCs. 

• The properties of the cell expressing the gene variant were studied using electrophysiology or Ca^"*" imaging. 

• The deviation between the variant cell property and the control cell property could be implemented in the applied 
neuron model [S2] as a change of a model parameter value. 

• The observed effect of the gene variant was not solely on the expression or ion channel density level. 

The last condition, which ruled out studies where the effect was only shown on channel density, was set due to the multitude of 
pathways that may contribute to such an effect [S3]. By contrast, the effects on channel kinetics (activation and inactivation 
threshold potentials, sensitivity to membrane potential, and opening/closing time constants) were expected to be more 
straightforwardly dependent on the way these channels are genetically coded. However, data on differences in ion channel 
expression bear an important aspect of the total ionic behavior as well, and should therefore be introduced in the future, 
although compensatory mechanisms may contribute to cancel the effects of changes in expression in vivo (cf. [S4]). 

Table 1 shows all studies that were found to obey the above conditions. An expansion of this table is given in Table SI, 
added with details about the effects of each variant and the underlying experimental data. The six first genes of Table 1 
belong to the calcium signaling toolkit. The genes CACNAlC, CACNB2, and CACNAlD, encode a or /3 subunits of HVA 
Ca^’*' channels, and CACNAlI gene encodes an a subunit of an LVA Ca^+ channel. Genes ATP2A2 and ATP2B2 encode 
intracellular Ca^’*' transporters, namely sarco/endoplasmic reticulum Ca^’*' ATPase (SERCA) and plasma membrane Ca^'*' 
ATPase (PMCA). The rest of the listed genes encode ion channels of other ionic species. SCNlA and SCN9A encode a 
subunits of Na"*" channels, which we considered to contribute to the currents I Nat and I Nap in the neuron model. KCNBl 
encodes a subunit that contributes to the slow, delayed rectifier K+ current in L5PCs [S5, S6], which is denoted by Ixp in the 
neuron model. KCNS3 encodes an electrically silent subunit, which yet together with the subunit encoded by KCNBl forms 
slowly inactivating K+ channels [S5, S7]: these channels contribute to the current Ikp as well. KCNN3 encodes a subunit 
for a Ca^^-activated small-conductance K+ channel (contributing to Isk), and HCNl encodes a subunit for a non-specific 
ion channel (contributing to It)- Not all of these genes have been shown to be expressed in thick-tufted L5PCs in specific, 
however, we confirmed in Allen Brain Atlas (mouse) and Human Protein Atlas that they all show expression in the cortex. 
Certain studies on CACNAlS fulfilled all other above conditions, but the gene is not expressed anywhere in the brain, and 
therefore we disregarded them. 

Some of the studies cited in Table 1 considered several variant types, and in these cases, the range of possible effects on 
model parameters is considered (see Table SI). If the range spanned both increasing and decreasing effects, the analyses were 
carried out for both end points of the range, otherwise only the maximal deviation from the control value was considered. 
Also, in case one of the end points of the range was very close to the control value, namely, at a distance less than 1 mV 
or less than 10% of the distance between control value and the other end point, only the endpoint with the larger deviation 
was considered. 

The applied model of an L5PC is relatively abundant in biological detail, yet it groups together many different types 
of currents under one model current. This represents a hindrance to carrying out detailed computational studies on the 
contribution of different genes to neuron firing behavior. As for our study, for this reason we could not make a distinction 
between the contributions of L, N-, P-, and Q-type channels (HVA Ca^+ current channels) to the activation of the SK current, 
which have been experimentally studied in e.g. [S8]. Taking into account this particular distinction would be an important 
extension to our modeling framework, as many of the SCZ-related Ca^"*" channel-encoding genes specifically encode L-type 
channel subunits {CACNAlC, CACNAlD, CACNAlS, CACNB2). Nevertheless, currently there are no L5PC models that 
would allow this. 

Scaling of gene variants 

Due to the polygenic nature of SCZ, it has been proposed that none of the SCZ-related SNPs can alone cause the disorder, 
and that the disease may only be induced when sufficiently many of them are represented. This paradigm was adopted 
into the computational study at hand as follows. We hrst identified the parameter changes in the L5PC model [S2] that 
corresponded to the effects of a certain genetic variant. We then ran simulations on the neuron model, exploring the effects 
of one variant separately on the neuron’s response to selected stimuli. If the variant altered the neural response dramatically, 
the genetic effect was scaled down (i.e. parameters were brought closer to the control values), so that there were no large 
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discrepancies between the control neuron and the mutant with the downscaled variant. The criteria for the downscaling were 
that no variant should alone change the number of action potentials to chosen stimuli or radically change the steady-state 
firing behavior of the neuron (control neuron referring to the model neuron exactly as described in [S2]). This approach was 
used in order to simulate the SCZ-related SNP effects, many of which are known to be subtle (cf. [S9]) in contrast to the 
large phenotypic effects caused by the variants of Table 1. 

To scale down the effects of the variants of Table SI, the following procedure was applied. Firstly, a “standard spiking 
behavior” was defined as a firing behavior that fulfils the following five conditions: 

(I) Exactly 4 spikes should be induced as a response to somatic square-current injection of 0.696 nA x 150 ms, 

(II) Exactly 1 spike should be induced as a response to a distal (620 pm from soma) alpha-shaped synaptic conductance 
(time constant 5 ms, max. amplitude 0.0612 pS), 

(III) Exactly 2 spikes should be induced as a response to a combined stimulus of somatic square-current injection (1.137 
nA X 10 ms) and distal synaptic conductance (time constant 5 ms, max. amplitude 0.100 pS, applied 5ms after the somatic 
pulse), 

(IV) The integrated difference between the f-I curves of the considered neuron and the control neuron should not be more 
than 10% of the integral of the control neuron f-I curve, and 

(V) The membrane potential limit cycle should not be too different from the control neuron limit cycle (dcdlci, IC2) < 600, 
see the section below for details). 

The above-mentioned amplitudes 0.696 nA, 0.0612 pS, and 1.137 nA were chosen such that the control neuron most stably 
produces the named numbers of spikes with the default parameters — most stably in the sense that an equal change in 
current amplitude on logarithmic scale is required in order to produce one spike more or to produce one spike less. The 
conditions I-III restrict the magnitudes of short-time responses of the neuron, while the conditions IV-V concern continuous, 
steady-state firing. Secondly, the spiking behavior of cells expressing the genetic variations listed in Table 1 under the same 
stimuli is simulated. In case one or more of the conditions I-V are violated, the effect is scaled down, all parameters in 
proportion, to a fraction c < 1 of the original effect where the violation is for the hrst time observed. These threshold effect 
factors c are listed in Table 2 for each variant together with the corresponding parameter changes, and the six variants that 
were used in Figures 2-6 are highlighted. As we do not know how small a parameter change effect should be to correspond 
to a small SNP effect, we consider variants with different scalings where the threshold effect factor c is multiplied with 
another factor e < 1. In this work, we consider the scaling factor values e = ^ and and we also display the effects of the 
corresponding opposite variants e = — 5 and — ^. 

For those unsealed variants that did not violate the conditions I-V, we sought for the threshold effect up to twice the 
original effect, i.e. c < 2. If the variant still obeyed conditions I-V, we considered the original variant the e = | variant and 
applied other scalings with respect to this choice. 

The model parameters affected by the variants include quantities of various roles and dimensions (mV, mM, ms, etc.), 
which calls for a careful consideration how to scale them properly. We chose to perform this such that those parameters that 
may receive both negative and positive values were scaled linearly, while the parameters that receive only positive values 
were scaled on the logarithmic scale. In practice, this means that the differences in offset and reverse potentials (I4ftm»j 
Eih) between control and variant neuron were expressed as an additive term (±a; mV), and this term x was multiplied 
by a factor c G [0,1] in the downscaling procedure. By contrast, the differences in all the other model parameters (t4ffc*, 
I4io», T*, c», 7 *) between control and variant neuron were expressed as a multiplication (xa;), where the downscaling caused 
this coefficient x to be exponentiated by the same factor c. This procedure permits a continuum of parameter changes in 
the range c G [0,1] that is directly applicable to amplified (c > 1) parameter changes as well. As an example, consider the 
mid-point of activation of the slow inactivating K+ channel, I4ffm,Kp = ~11 mV. One of the variants of KCNBl gene (fifth 
in Table 2) shifted this value by -1-14 mV, which gives the new value Kp = ^ ^ case, it is not possible scale 

the parameter logarithmically, as the factor x = Kp/^offm,Kp to be exponentiated is negative. As a contrary example, 
the same variant changes the time constants of activation, Tm*,Kpj by —61%. In principle, this effect could be scaled both 
logarithmically and linearly in the range c G [0,1]. However, if the unsealed variant does not cause a violation of conditions 
I-V, we would study the effects of an “upscaled” variant all the way until c = 2. This would not be possible using linear 
scaling, as we would have to decrease the time constants of activation by 122 %, leading to negative time constants. 


A distance metric for membrane potential limit cycles 


A membrane potential limit cycle can be described as a 1-dimensional manifold in the space CC = V x dV, where V is the 
space of possible values for membrane potentials and dV is the space of possible values for time derivatives of membrane 
potential. Due to the difference in units for x and y axis, there is no obvious metric for this space, but such can easily be 
constructed. We define a bijection / : £C —)■ as 


fiV,dV) 


V 


dV 


1.0 mV’ 7.026 mV/ms 


where the constant 7.026 was chosen such that the control neuron limit cycle spans an equal range on x and y axes when 
the soma is given a DC of amplitude 1.0 nA. Now, a distance metric in d : CC x CC can be dehned as 


d{x,y) = ||/(x) - f{y)\\ 2 , 


where jj • II 2 is a Euclidean norm. Furthermore, the distance of a point x G CC from a limit cycle Ic C CC can be dehned as 


ddxjc) = min||/(a;) - f{y)\\. 

y£lc 
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To evaluate the difference between two limit cycles ^ci, lc 2 C CC^ we define 


dcc(^Ci, lc2^ 


Ici 


dc{x, lc2)dx - 


lC2 


dc{x, lci)dx. 


Using the average of the two integrals makes sure that no part of either limit cycle is ignored: If only integral over points 
of Ici was used, the limit cycle lc 2 could have any “extra” part, for example an unusually long and slow hyperpolarization 
period or minor suboscillations, that might largely be disregarded in the integration as there are always some other points 
in /c 2 that lie nearer to Ici. 
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Table SI: Table of genetic variants and their maximal effects on our neuron model parameters. The first column 
of the table shows the gene whose variant is studied. The second column shows the model parameters that are affected by 
the variant, as the third column shows the size of the effect. The fourth column gives the PubMed index of the article, 
where the effects of the variant are measured. The hfth column names the type of variant used, while the sixth and seventh 
columns show the cell type in which the effects are measured and the animals species used in the study. The final column 
may give extra information. The variants are listed in the order of genes that they affect. Asterisks represent an identical 
change to a set of variables as follows (applicable to those variables that exist in the considered ion channel type): offm* 
refers to variables VoSm, Voffma and VoBmh', offh* refers to variables VoSh, fbffha and Vbffhb; slom* refers to variables f4ioma 
and Psiomb; sloh* refers to variables f4ioha and Psiohb; taum* refers to variables r^a, Tmb, Tmmin, Tmdiff, Tmdiffi and rmdiff 2 ; 
tauh* refers to variables Tha, Thb, Thmin, Thdiff, Thmean, Thdiffi and rhdiff 2 - 


Gene 

Model parameter 

Effect 

Paper 

Type of variant 

Cell type 

Animal 

Notes 

CACNAIC 

of fm*_CciHVA 
offh*_CaHVA 

-25.9..-1.4 mV 

-27.0..-3.8 mV 

pubmed/19265197 

L429T, L434T, S435T, S435A 
S435P 

TSA201 

human 

-These represent variant of the 
IS6 segment 

CACNAIC 

offm*_CaHVA 

offh*_CaHVA 

-37.3..-9.7 mV 
-30.0..-11.8 mV 

pubmed/19265197 

L779T, I781T, I781P 

TSA201 

human 

-These are variants of IIS6 
segment. Double (IIS6+IS6) mutant 
seems to have additive effect 

CACNB2 

offh*.CaHVA 

sloh*_CciHVA 

-5.2 mV 
-317, 

pubmed/19358333 

Tlll-mutant 

TSA201 

humcin 

-Small changes shown in offma, 
offmb and tau ignored 

CACNB2 

taum*_CaHVA 

+707, 

pubmed/7723731 

A1B2 vs A1 alone 

HEK293 

human/raouse 


CACNB2 

offm*_CaHVA 

offh*_CaHVA 

taum*_CciHVA 

tauh*_CciHVA 

-4.9..+4.9raV 

-5.1..+5.IraV 

-407,. .+687. 

-407,. .+667. 

pubmed/19723630 

N1 vs N3 vs N4 vs N5 

HEK293 

humcin/raouse 

-Changes to voltage-dependence 
of taus ignored. 

CACNAID 

offm*_CaHVA 

slom*_CaHVA 

offh*_CaHVA 

sloh*_CaHVA 

tauh*_CaHVA 

-10.9..-8.5mV 
-27. .-137, 

-3.0..+3.5raV 
-12. .-197, 

+257, 

pubmed/21998309 

pubmed/21998310 

42A splices trainsfected 
(Alternative splicing at C- 
terminus) 

TSA201/HEK293 

human 

-The variants are expressed in 
human and mouse brains; hence, 
the effect of genetic varicinces 
could be of opposite sign with 
respect to control 

CACNAID 

of fm*_CciHVA 
slom*_CciHVA 
offh*_CaHVA 
sloh*_CaHVA 
tauh*_CaHVA 

-10.6..+3.4mV 
-20. .+127, 

-5.3..+1.2raV 

-34. .-87, 

-287, 

pubmed/21998309 
pubmed/21998310 

43S splices trzinsfected 
(Alternative splicing at C- 
terminus) 

TSA201/HEK293 

humcin 

-The variants are expressed in 
humcin and mouse brains; hence, 
the effect of genetic variances 
could be of opposite sign with 
respect to control 

CACNAID 

offm*_CaHVA 

slom*_CciHVA 

tauh*_CaHVA 

+3.5..+6.6raV 

-25. .+197, 

-50. .+127, 

pubmed/20951705 

pubmed/21054386 

CaVl.3 (-/-) vs control 

AV node cells/ 
chromaffin 

cells 

mouse 

-Some time constcints compared 
between single tau fits (as 
double tau fits not always well 
fitted) 

CACNAII 

offma.CaLVA 
offha.CaLVA 
taum*_CaLVA 
tauh*_CaLVA 

-0.2..+1.3raV 

-0.5..+1.6raV 
-13. .+457, 

-20. . +87, 

pubmed/15254077 

Alternative splicing of exons 

9 and 33 

HEK293 

human 

-Maximum of effects at -40 or 

0 mV on kinetics considered 
-Chcinges in slope minuscule 

ATP2A2 

gamma.CciDyn 

-30. . -407, 

pubmed/10970890 

Heter. null mutation 

myocytes 

mouse 

-Compensating mechanisms exist 
that prevent larger effects 

ATP2B2 

decay.CciDyn 

+ 15. .+1137, 

pubmed/22789621 

Heter. knockout 

Purkinje cells 

mouse 

-Compensatory mechanisms exist 
that prevent larger effects 

ATP2B2 

decay.CaDyn 

minCai.CaDyn 

+32. .+507, 

+407, 

pubmed/21232211 

Homozyg. knockout 

Purkinje cells 

mouse 


ATP2B2 

decay.CciDyn 

+53. .+3457, 

pubmed/17234811 

G283S-,G293S-rautant 

CHO cells 

human/raouse 

-+537, for G293S, +3457, for G283S 

SCNIA 

offm_Nat 
offh_Nat 
slom.Nat 
sloh_Nat 

-0.3mV 
+5.OmV 

+ 157, 

+237, 

pubmed/18632931 

FHM mutation Q1489K 

Neocortex 

(cultured) 

human/rat 

-Slow inactivation could not be 

studied in detail in neurons 

SCNIA 

offm_Nat 
offh_Nat 
slom.Nat 
sloh_Nat 

+2.8mV 
+9.6mV 

-1.67. 

+4.27. 

pubmed/17397047 

FHM mutation L1649Q 

TSA201 

human 

-Electrophysiology done with the 
corresponding mutation L1636Q in 
the homologous SCN5A gene due to 
instabilities in recomb, bacteria 

SCN9A 

offh*_Nap 

+6.8mV 

pubmed/22136189 

I228M NaVl.7 variant 

HEK293 

human 

-Differences in activation and 
fast inactivation not significant 

SCN9A 

offh*_Nap 
sloh_Nap 
offm_Nat 
offh_Nat 
sloh_Nat 

+3.5mV 

-457, 

-7.ImV 

+17.0raV 

-317, 

pubmed/18945915 

A1632E NaVl.7 mutation 

HEK293 

human 

-Slow inactivation interpreted as 
inactivation of persistent (Nap) 
current and fast inactivation as 
inactivation of transient (Nat) 
current 

SCN9A 

offm_Nat 
offh_Nat 

-9.ImV 

+3.ImV 

pubmed/16392115 

L858F NaVl.7 mutation 

HEK293 

human 

-Activation and inact. curves 
closer to Nat them Nap currents 

SCN9A 

offm.Nat 
offh_Nat 

-7.6mV 
+4.3mV 

pubmed/15958509 

F1449V NaVl.7 mutation 

HEK293 

human 

-Activation and inact. curves 
closer to Nat them Nap currents 

KCNS3 

taum*_Kp 

tauh*_Kp 

sloh.Kp 

0. .+1007, 

+50. .+1507. 

-507, 

pubmed/10484328 

hKv2.l-(G4S)3-hKv9 fusion 
inserted (or hKv9 inserted 
explicitly) 

HEK293 

human 

-Fits not carried out in the 
paper, approximate values used. 

KCNBl 

offm_Kp 

offh_Kp 

slom.Kp 

sloh_Kp 

taum*_Kp 

tauh*_Kp 

+5mV 

+3mV 
+ 117, 

-147, 

-507, 

-477, 

pubmed/21455829 

T203K mutcint 

HEK293 

human/raouse 

-Double mutations studied as 
well, but they are not included 
here 

KCNBl 

of fm_Kp 
offh_Kp 
slom.Kp 
sloh_Kp 
taum*_Kp 

+ lmV 

-6mV 

+227, 

+07, 

-117, 

pubmed/21455829 

T203D mutcint 

HEK293 

humcin/raouse 


























tauh*_Kp "13% 


KCNBl 

offm_Kp 

offh_Kp 

slom.Kp 

sloh_Kp 

taum*_Kp 

tauh*_Kp 

+6mV 

-8mV 

+33% 

+0% 

-50% 

-13% 


pubmed/21455829 

S347K mutcint 

HEK293 

humcin/mouse 


KCNBl 

offm_Kp 

offh_Kp 

slom.Kp 

sloh_Kp 

taum*_Kp 

tauh*_Kp 

-28mV 
-27mV 
+ 11% 
-29% 

+ 13% 

+ 127% 


pubmed/21455829 

S347D mutcint 

HEK293 

human/raouse 


KCNBl 

offm_Kp 

offh_Kp 

slom.Kp 

sloh_Kp 

taum*_Kp 

tauh*_Kp 

+ 14mV 

-21mV 

+ 100% 

+0% 

-61% 

+20% 


pubmed/21455829 

T203W mutcint 

HEK293 

human/raouse 


KCNBl 

offm_Kp 

offh_Kp 

slom.Kp 

sloh.Kp 

taum*_Kp 

tauh*_Kp 

-13mV 

-13mV 

+33% 

-29% 

-5% 

+413% 


pubmed/21455829 

S347W mutcint 

HEK293 

human/raouse 


KCNN3 

offc.SK 

sloc.SK 

-14% 

+24% 


pubmed/14978258 

hSK3_ex4 isoform 

TSA 

human 

-Activation curve constructed a 
bit differently from Hay model 

HCNl 

offm*_Ih 
slom*_Ih 

-2.1.. 
-12. 

-26.5mV 

•36% 

pubmed/17185333 

D135W, D135H, D135N mutants 

HEK293 

humcin/raouse 

-Other mutations studied as well, 
but they change the Ih current 
too dramatically 
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Supplemental results 

Here, we extended the analyses carried out for the variants of Figures 2-6 by considering additional variants. Figures SI 
and S2 show a compilation of statistics related to Figures 3-6 for all variants of Table 2: Figure SI illustrates the range of 
steady-state Ca^’*' concentration oscillation, while Figure S2 shows the mean firing rate across the DC amplitudes used in 
the f-I curves, the temporal sensitivities of the Ca^“'" spike during up and down states, and the time window for an efficient 
suppression of a second apical stimulus. In Figure SI, the results are additionally plotted for the unsealed variants for 
reference, while in Figure S2, only unsealed variants are considered. The results of Figure S2 for the positively scaled (e = |) 
variant are also summarized in Table S2. Except for genes KCNS3 and KCNBl, variants of all genes showed a notable 
impact on some of these statistics. 

Changes in intracellular [Ca^+J may be functionally relevant for many neural processes. Intracellular [Ca^+J plays an 
important role in the activation of various second messenger cascades, which control a variety of cellular processes [S4, 
S10-S13], but these aspects are out of the scope of the current work. A more immediate effect of intracellular [Ca^+J is the 
activation of Ca^+-dependent K+ -currents, such as the small-conductance K+ current (Isk)- This current was included 
in the L5PC-model, and it is affected by the KCNN gene family. In Figure 2, the effects of the KCNN3 variant were 
undetectable. This was because Isk currents are large only when both intracellular [Ca^’*'] and the membrane potential are 
elevated — although the currents only need high [Ca^'*'] to be activated, they can be large only when the membrane potential 
is far enough from the K'*' reversal potential, which is relatively near to the membrane resting potential. The effects of the 
KCNN3 variant could be detected when a more enduring stimulus than a brief somatic pulse was applied, as done in Figures 
3-6: these effects can also be observed in Figures SI and S2. 

Most of the variants that gave positive shift in [Ca^+](Figure SI), also resulted in reduced firing rates (Figure 3). There 
was a correlation coefficient of -0.70 between the maximal [Ca^'*'] value in the limit cycle and the integral of the measured f-I 
curve. An exception to this trend was the KCNN3 variant, where the threshold [Ca^+J was lower (in the positively scaled 
variants) than in the control neuron, see Table 2. When stimulated with the DC of amplitude 1.2 nA (Figure 4), this lower 
threshold caused the SK current to be upheld for a longer time after the spike than in the control neuron. This caused a 
delay in the timing of the next spike, allowing the [Ca^+] to drop to a lower value before the next spike, and hence the 
oscillation of the [Ca^^] remained at an overall lower level than in the control neuron (data not shown). 

Combinations of variants may cause large alterations of neuron excitability 

The effects of the variants on neuron excitability, as explored in Figures 2-6, were relatively small when the variants were 
studied in isolation (this was also the intention behind the downscaling process). However, combinations of variants of 
different genes may have additive effects on the neuron excitability. We investigated the effect that combinations of several 
of the variants could have. Figure S3 shows a selected example of how a combination of variants altered the neuron firing 
properties considered in Figures 2-6. Although this version of the model neuron consists of downscaled variants, the properties 
of neuron excitability were remarkably modified. 

In Figure S3, the selection was based on the property of inhibiting the second apical stimulus as shown in Figure 6: We 
combined such downscaled variants that had a lower threshold curve than the control neuron, i.e., variants that made it 
easier for the second apical stimulus to induce a spike, however picking maximum one variant per gene. See Table 2 for 
this selection. Large impacts of the combination of variants could be seen in both Ca^’*' responses (Figure S3A-B, S3D), 
temporal profiles of Ca^^ spike generation (Figure S3E), and naturally the inhibition of a second stimulus (Figure S3F). 
The change in gain (Figure S3E) was nevertheless minor: Although all variants of the combination caused an increase in 
gain when in isolation (except for the CACNAII variant which had a minuscule effect and KCNN3 variant which had an 
ambiguous effect), these increases in gain were in most cases caused by an overall decrease in Ca^’*' currents. These changes 
were compensated for by a decrease of threshold [Ca^+] of the SK current that is characteristic of the KCNN3 variant, and 
hence the resulting f-I curve was relatively close to that of the control neuron. 

The effect of the combination of variants on the Ca^+ spike generation in Figure S3E is notable. The sensitivity of the 
Ca^’*' spike magnitude to the ISI could serve as a coincidence detection mechanism of synaptic inputs projected to different 
parts of the L5PC and thus form an important element in cortical information processing. A widening of the region of 
ISIs that produce a large Ca^+ spike could weaken the coincidence detection ability of the neuron [S14], while an extensive 
narrowing of the curve would limit the regime for [Ca^“'"] spike generation [S15] and hence contribute with a lesser total 
excitability. In the up-state protocol, the combination of e = j variants only produce minor Cs?'^ spikes, and their temporal 
ISI window is narrowed down, while the combination of e = ^ variants does not produce a Ca^"*" spike at all with the 
considered stimuli. 

Due to lack of data on the interplay between different variants of a single gene, we only implemented a maximum of 
one variant per gene when we simulated combinations of variants. There is, however, evidence that SNPs or other variants 
located in different parts of the same gene may have additive effects on the ion channel kinetics, as shown in the case of 
mutations in different parts of the CACNAIC gene [S16]. Hence, given detailed information on such effects in other genes 
as well, the predictive power of our method could still be improved by allowing more diverse combinations of the same set 
of variants. 


Robustness of the analyses 

To study the robustness of our findings, we performed similar analyses in a neuron model using an alternative set of fitted 
model parameters. We implemented the alternative neuron model, “Model 1”, as given in the first column of Table SI in [S2]. 
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The scaling of the variants was redone in this neuron model, and the variant neuron behaviour was quantified in a similar 
manner as done for the primary fit. Figure S4 shows the compilation of results obtained from the alternative model. The 
results are qualitatively similar for the majority of the variants. The only observable exceptions are the 12th variant in the 
third group of CACNB2 variants and the second variant in the second group of CACNAlD variants, which caused a large 
depression of the Ca^^ spike in the up-state protocol of the primary neuron model, shown as a missing bar in Figure S2C. 
By contrast, the dendritic membrane potential response in the same variants of the alternative neuron model did exceed the 
set threshold, leading to an observed (but narrowed-down) temporal window in Figure S4C. 

We further analyzed the robustness of our findings using the primary neuron model fit but in an alternative neuron 
morphology. We applied and downscaled our variants using the cell ^2 morphology of [S2], and performed the same analyses 
as done with the primary neuron morphology, cell ^1. The firing rates were in general smaller in cell ^2: A somatic DC of 
amplitude 1.2 nA caused the control neuron to fire with a frequency of 15.4 Hz, whereas the corresponding rate in cell ^1 
was 16.9 Hz. The up- and down-state protocols did not produce responses directly comparable those of cell #1, and hence, 
we omitted these analyses in cell #2. By contrast, the property of inhibiting a second apical stimulus was preserved in cell 
^2, although in a slightly different regime. In cell #2, a nearly doubled maximal conductance w.r.t. that in cell #1 was 
required to make the neuron spike when the synaptic inputs were uniformly distributed along the apical dendrite {gth = 5.0 
pS vs. 2.8 pS in cell #1). Accordingly, the ratio Cg between the threshold conductance for eliciting a second spike and that 
for eliciting the first spike {gth) was lower: In the control case of cell #2, the factor Cg varied from 0.44 to 1.91, while the 
corresponding range in cell #1 was from 0.51 to 5.13. To compensate for this effect, we analyzed the PPI temporal window 
in terms of a lower threshold factor, namely, we measured the first and last ISI for which Cg exceeded the value 1.5 (3.0 in 
cell #1). 

The compilation of results obtained from simulations with cell ^2 are shown in Figure S5. The amplitudes of variant 
effects are different than in cell #1 (see Figure S2), but the direction of the variant effects are the same as in cell #1. 
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Figure SI: Ranges of Ca^+ concentration in simulations of steady-state firing in cells with different genetic 
variants implemented. The blue bars (and the light blue horizontal lines) show the minimum and maximum Ca^'*' 
concentration as a response to a 1.2 nA DC injection in a control neuron, while the other colors show the corresponding 
quantities in neurons with the variants from Table 2 implemented. The red bars in the upper plot show the effects of the 
unsealed variants, while the bars in the middle plot show the effects of the e = | (magenta) and e = —^ (green) downscaled 
variants. The bars in the lower plot show the effect of minimally scaled (dark gray: e = j, cyan: e = — j) variants. 
The missing red bars among CACNAlC and SCN9A variants represent cases where the unsealed variant expresses large 
depolarization even without an input, and the membrane potential repolarization is missing or too weak to maintain the 
spiking behavior. The variants are based on previous experimental studies [S7, S16-S37] and they are ordered as in Table 2. 
The variants marked with f were used in Figures 2-6. 
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Figure S2: Various statistics from neurons with different downscaled genetic variants implemented. A: The 

average steady-state spiking frequency (see Figure 3) when a somatic DC of amplitude 0.35-1.4 nA was injected. Data 
plotted for each variant in Table 2. B: The down-state temporal window span (see Figure 5). The plotted range corresponds 
to the range of ISIs for which the maximum membrane potential at the apical dendrite was larger than -14.05 mV (the 
midpoint between the highest and lowest apical dendrite membrane potentials across all ISIs for the control neuron). C: The 
up-state temporal window span (see Figure 5) of ISIs that produced a peak apical dendrite membrane potential larger than 
-23.21 mV, illustrated in a similar manner as in panel B. D; The temporal window for large single-cell prepulse inhibition 
(PPI, see Figure 6 ). The plotted ranges correspond to the ISIs for which a larger than three-fold second apical stimulus 
was required to induce an extra spike. In all panels, the blue bars (and the light blue horizontal lines) correspond to the 
control neuron values, while the magenta and green bars correspond to the e = 5 (magenta) and e = — 5 (green) downscaled 
variants, and the dark gray and cyan bars correspond to the e = j (dark gray) and e = — | (cyan) downscaled variants. The 
missing bars in panel B among SCN9A variants correspond to cases where the somatic pulse alone was sufficient to produce a 
large Ca^'*' spike (apical dendrite membrane potential larger than -14.05 mV). The missing bars in panel C among CACNB2 
and CACNAID variants correspond to cases where no ISI produced a large Ca^+ spike (i.e., the apical dendrite membrane 
potential always remained smaller than -23.21 mV). The variants marked with f were used in Figures 2-6. 
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Figure S3: Combined genetic variants can radically change the neuron response. A combination of six downscaled 
variants from Table 2 {CACNB2 [S19], CACNAlD [S22, S23], CACNAlI [S24], ATP2A2 [S25], and KCNN3 [S36]) was 
applied. These variants were chosen as an example combination as they all showed weakened prepulse inhibition (see Figure 
6 ). The gray curves show the properties of a neuron with the combination of e = j downscaled variants, and the magenta 
curves represent the combination of e = ^ downscaled variants. A-B: [Ca^+] response to short somatic stimulus, shown 
in time series (A) and phase plane (B). See Figure 2B and 2C for details. C: f-I curves. See Figure 3 for details. D: 
[Ca^+] response to prolonged stimulus. See Figure 4 for details. E: Integration of apical and somatic stimulus during up 
and down states. See Figure 5 for details. The combination of e = ^ variants did not respond with a strong Ca^+ spike for 
any inter-stimulus interval. F: Threshold conductance factor for a second spike. See Figure 6 for details. A strong second 
stimulus with four-fold total conductance with respect to the threshold conductance would always make the neuron with 
combined e = | variants spike, while the neuron with combined e = j variants would inhibit it for ISIS[49 ms, 61 ms], and 
the control neuron would inhibit it for ISIS [49 ms, 71 msj. 
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Table S2: Effects of positively scaled variants on steady-state firing, integration of somatic and apical inpnts, 
and inhibition of a second apical stimulus. The shown effects are extracted from the properties of the e = | variant 
in Figure S2. The variants are ordered as in Table 2, and the hrst column shows the name of the underlying gene. The 
second column shows whether the average firing rate (see Figure S2A) of the variant is increased (’+’) or decreased 
with respect to the control neuron. A “minuscule” effect refers to cases where the effect was smaller than 0.5% of the control 
neuron average firing rate. The third and fourth columns show how the down and up-state windows (see Figure S2B-C) are 
affected in the variant. A “broader” window refers to an effect where the variant window starts at an ISI smaller than the 
start of the control window and ends at an ISI greater than the end of the control window. By contrast, in a “narrower” 
case the window starts at an ISI greater than the start of the control window and ends at a smaller ISI than the end of the 
control window. If both the starting point and end point of the window are moved to the same direction with respect to the 
control window, the table entry is “moved forward” (toward smaller ISIs) or “delayed” (toward larger ISIs). Cases where 
both end points moved less than 5% of the width of the temporal window in a control neuron are labeled as “minuscule” 
effects, regardless the direction to which the end points moved. Finally, the value “N/A” denotes cases where the temporal 
window could not be defined. The fifth column shows the corresponding properties for the PPI window (see Figure S2D). 
The variants marked with f are the ones that were used in Figures 2-6. The horizontal lines separate variants that are based 
on different sets of experimental data, see Table 1. 
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Figure S4: Various statistics from neurons with alternative set of fitted model parameters and different 
downscaled genetic variants implemented. The maximal conductances and parameters controlling Ca^'*' dynamics 
from Table SI in [S2] (“Model 1”) were applied, while other parameters were kept fixed. A: The average steady-state spiking 
frequency when a somatic DC of amplitude 0.35-1.4 uA was injected. See Figure S2A for the corresponding results in 
the neuron model with default parameters. B: The down-state temporal window span. See Figure S2B. C: The up-state 
temporal window span. See Figure S2C. D: The temporal window for large single-cell prepulse inhibition. See Figure S2D. 
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Figure S5: Average firing frequency and PPI statistics from neurons with alternative morphology and different 
downscaled genetic variants implemented. The cell morphology #2 of [S2] was used with the default model parameters. 
A: The average steady-state spiking frequency when a somatic DC of amplitude 0.35-1.4 nA was injected. See Figure S2A 
for the corresponding results in cell #1. B: The temporal window for large single-cell prepulse inhibition. See Figure S2D. 
The plotted ranges correspond to the ISIs for which a larger than 1.5-fold second apical stimulus was required to induce an 
extra spike. 
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